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Abstract  A  method  is  presented  for  the  modeling  of 
brittle  elastic  fracture  which  combines  peridynamics 
and  a  finite  difference  method  to  mitigate  the  wave 
dispersion  properties  of  peridynamics.  Essentially,  a 
finite  difference  method  is  used  in  the  bulk  for  wave 
propagation  modeling,  while  peridynamics  is  automat¬ 
ically  inserted  in  high  strain  areas  to  model  crack 
initiation  and  growth.  The  dispersion  properties  of 
finite  difference  methods  and  discretized  peridynam¬ 
ics  are  reviewed  and  the  interface  refiection  properties 
between  the  two  regions  are  investigated.  Results  show 
that  the  augmented  method  can  improve  the  modeling 
of  wave  propagation  and  boundary  conditions.  In  addi¬ 
tion,  the  numerical  stress  intensity  factor  computed  at  a 
crack  tip  shows  reduced  oscillations  in  the  augmented 
method,  likely  due  to  the  improved  dispersion  proper¬ 
ties  of  the  bulk.  Dynamic  fracture  simulations  show  a 
difference  in  crack  paths  between  the  methods. 

Keywords  Einite  difference  •  Peridynamics  • 
Eracture  •  Wave  propagation 

1  Introduction 

Computational  mechanics  problems  involving  the  frac¬ 
ture  and  failure  of  materials  are  difficult  to  solve  using 
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partial  differential  equation  formulations  as  any  dis¬ 
continuity  introduced  by  such  failure  causes  difficul¬ 
ties  in  both  the  computation  and  formulation.  Peridy¬ 
namics  proposes  to  address  this  issue  by  replacing  the 
spatial  differential  operators  in  the  standard  continuum 
mechanics  equations  with  integral  operators  that  con¬ 
verge  to  those  differential  operators  under  certain  con¬ 
ditions  (Silling  2000;  Silling  et  al.  2007;  Emmrich  and 
Weckner  2006).  When  discretized  and  coupled  with 
a  damage  model,  peridynamics  has  been  shown  to  be 
capable  of  matching  a  number  of  experimental  results 
in  both  brittle  and  ductile  fracture  problems  (Silling 
and  Askari  2005;  Ha  and  Bobaru  2010,  2011;  Wild¬ 
man  and  Gazonas  2013;  Poster  et  al.  2010).  In  addi¬ 
tion,  peridynamics  has  been  applied  to  other  fields  of 
physics,  such  as  heat  transfer  (Bobaru  and  Duangpa- 
nya  2010)  and  electromigration  (Gerstle  et  al.  2008). 
While  peridynamics  has  seen  success  in  matching  frac¬ 
ture  patterns  observed  in  experiments,  issues  remain, 
such  as  its  behavior  at  material  boundaries  and  its  dis¬ 
persion  characteristics  (Weckner  and  Abeyaratne  2005 ; 
Mikata  2012).  Here,  we  attempt  to  address  dispersion 
in  peridynamic  simulations  by  coupling  with  a  finite 
difference  solver  to  compute  internal  forces. 

Einite  difference  methods  have  a  long  history,  begin¬ 
ning  with  a  proof  of  their  stability  properties  for  hyper¬ 
bolic  partial  differential  equations  (Courant  et  al.  1928) 
and  similar  methods  exist  in  many  fields;  a  short 
list  includes  electromagnetics  (Alterman  and  Karal 
1968;  Zahradnik  and  Hron  1992),  elasticity  (Yee  1966; 
Tafiove  andBrodwin  1975),  and  diffusion  (Brian  1961). 
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Their  first  implementations  relied  on  rectangular,  uni¬ 
form  grids,  though  finite  difference  methods  have  been 
developed  for  unstructured  grids  and  curved  regions 
(Hesthaven  and  Warburton  2002). 

A  finite  difference  method  (FDM)  will  be  aug¬ 
mented  with  peridynamics  in  this  work,  as  the  most 
basic  implementations  of  FD  methods  are  point-based 
on  uniform  grids,  much  like  standard  peridynamic  dis¬ 
cretizations.  The  compatibility  of  the  spatial  discretiza¬ 
tion  allows  for  straightforward  switching  between 
the  two  methodologies.  Essentially,  all  that  changes 
between  the  two  formulations  is  the  calculation  of  inter¬ 
nal  forces:  in  FD  methods,  an  approximation  of  the 
divergence  of  the  stress  tensor  is  used,  and  in  peridy¬ 
namics,  the  integral  of  a  micro-force  function  is  used 
over  some  finite  horizon.  In  fact,  peridynamics  can  be 
viewed  as  a  sum  of  continuous  differences,  and  when 
discretized  appears  much  like  a  weighted  sum  of  finite 
differences.  In  addition,  peridynamics  is  known  to  be 
less  accurate  on  boundaries  so  the  finite  difference  dis¬ 
cretization  can  be  used  to  implement  boundary  condi¬ 
tions  in  undamaged  regions. 

Adaptive  methods  for  peridynamics  have  been  intro¬ 
duced  to  improve  its  accuracy  (Bobaru  et  al.  2009; 
Bobaru  and  Ha  2011).  These  methods  use  an  adap¬ 
tive  grid  to  refine  the  solution  around  crack  tips  and 
other  stress  concentrations.  It  has  been  shown  that  peri¬ 
dynamics  is  more  accurate  and  converges  to  the  con¬ 
tinuum  mechanics  formulation  as  the  horizon  shrinks. 
Of  course,  a  shrinking  horizon  must  be  coupled  with 
a  refined  grid,  hence  the  adaptive  refinement  method. 
Here,  we  attempt  to  gain  more  accurate  wave  propa¬ 
gation  using  a  different  formulation  in  the  bulk,  rather 
than  adaptive  refinement.  As  will  be  shown,  a  more 
accurate  dispersion  relation  leads  to  more  accurate 
stress  concentrations  around  crack  tips  than  peridynam¬ 
ics  alone,  without  the  need  for  adaptive  refinement. 

Currently,  hybrid  peridynamic  s/finite  element  meth¬ 
ods  have  been  developed  (Lubineau  et  al.  2012;  Azdoud 
et  al.  2013,  2014).  The  combination  of  peridynam¬ 
ics  with  a  local  finite  element  method  has  the  same 
goal  as  this  work,  though  current  literature  is  restricted 
to  quasi-static  problems  and  the  method  requires  the 
local  and  non-local  regions  be  set  a  priori.  In  addi¬ 
tion,  Seleson  et  al.  (2013)  describe  a  hybrid  finite  dif- 
ference/peridynamics  model  using  a  blending  function 
in  ID,  though  again,  the  focus  is  on  preset,  unchang¬ 
ing  local/nonlocal  regions.  In  contrast,  this  work  will 
focus  on  the  dynamic  insertion  of  peridynamic  nodes 


into  a  finite  difference  domain.  Previous  approaches 
have  used  blending  functions  to  ease  the  reflections 
from  these  disparate  computational  domains,  though 
here  the  focus  will  be  on  dynamic  insertion. 

The  remainder  of  the  paper  proceeds  as  follows: 
In  Sect.  2,  peridynamics  is  summarized  and  a  specific 
finite  difference  method  is  formulated.  In  addition, 
the  dispersion  properties  of  both  discretizations  are 
reviewed,  and  finally  a  hybrid  approach  is  detailed.  In 
Sect.  3,  numerical  results  are  given,  including  a  demon¬ 
stration  of  the  dispersion  properties  of  both  methods, 
interface  refiections  between  the  two  computational 
regions,  a  comparison  of  stress  intensity  factors  at  a 
stationary  crack  tip,  and  finally  a  dynamic  crack  prop¬ 
agation  problem.  Section  4  concludes  the  paper. 

2  Formulation 

In  this  section,  first  peridynamics  and  FDM  are  sum¬ 
marized,  followed  by  a  discussion  of  their  numerical 
dispersion  properties.  Subsequently,  a  hybrid  method 
is  presented,  which  essentially  uses  finite  difference 
approximations  on  discretization  nodes  undergoing 
small  strain,  and  peridynamics  in  areas  of  high 
strain. 

2.1  Peridynamics 

To  facilitate  comparisons  of  dispersion  characteristics, 
the  peridynamics  formulation  is  the  linearization  of 
the  standard  bond-based  method  described  in  (Silling 
2000).  Bond-based  peridynamics  can  be  stated  as  the 
integro-differential  relation 

92  r 

P — yU  =  /  f  (u'  —  U,  —  x)  +  b,  (1) 

where  bold  type  indicates  a  vector,  u  =  u(x,t),u'  = 
,  p  is  the  density  (assumed  to  be  homoge¬ 
neous),  Hx  is  the  horizon  or  region  of  influence  at  point 
X  and  is  typically  a  sphere  of  radius  5 ,  b  is  a  body  force, 
and  f  is  a  peridynamic  force  function.  The  linearized 
version  of  Eq.  (1)  can  then  be  stated  as 

p^u=  /  c^vdx'  +  h,  (2) 

9?^  Jn, 

where|  =  x'— x,  t)  =  u'—u  and  c  is  a  material  constant 
that  can  be  related  to  classical  elasticity  constitutive 
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parameters.  For  example,  in  3D,  c  is  related  to  Young’s 
modulus  E  and  Poisson’s  ratio  v  as 
6E 

^  ^  7rS^{l-2v) 

where  v  =  1/4,  in  2D  plane  stress  we  have 


6E 

7t8^  (1  —  v)  ’ 


(4) 


where  v  =  1/3  and  in  2D  plane  strain  we  have 


6A. 

7t8^V  ’ 


(5) 


where  v  =  i/4  and  A.  is  Lame’s  first  parameter,  all  of 
which  are  derived  via  matching  the  elastic  strain  energy 
of  the  classical  and  peridynamic  formulations  (Ha  and 
Bobaru  2010). 

A  standard,  point-based  discretization  of  Eq.  (2),  is 
used  and  is  given  by 


Nm 

Z' 

n=l 


^mn  ^ 

nrr 

\Smn  I 


^mn  H” 


(6) 


where  Nm  is  the  number  of  nodes  in  the  horizon  Hx  of 
node  Xm,  ^mn  t)  -  U  (x^,  t), 

and  Vmn  is  the  volume  (or  area  in  2D)  of  the  node 
associated  with  node  m  (Silling  and  Askari  2005).  This 
distinction  on  the  volume  is  necessary  because  exact 
integration  weights  that  account  for  the  edge  of  the 
horizon  are  used.  In  other  words,  the  regions  or  cells 
associated  with  nodes  in  the  reference  configuration  are 
assumed  to  be  rectangular  (with  sides  given  by  the  dis¬ 
cretization  sizes  Ax  and  Ay),  and  if  the  entire  rectangle 
resides  within  a  node’s  horizon  then  its  area  is  simply 
the  area  of  the  square  (in  2D).  If  the  region  intersects 
the  horizon,  then  its  partial  area  is  computed  exactly 
by  using  the  divergence  theorem  to  convert  the  area 
integral  into  a  line  integral  that  encloses  the  integration 
region.  The  line  integral  will  then  be  a  collection  of 
zero  or  more  straight  line  segments  and  a  circular  arc, 
all  of  which  can  be  computed  exactly. 

The  temporal  discretization  is  a  velocity  Verlet 
method,  stated  as 

y^+l/2  ^yk 

u^+1 

_  (F  +  b)  (7) 

^  , 

P 

yk+i  ^  ^+1/2 


where  a  =  9^u/a^2  acceleration,  v  =  ^^/dt  is  the 
velocity,  F  is  the  internal  force  given  as  the  first  term 
of  the  right  hand  side  of  Eq.  (6),  b  is  a  given  body  force 
(that  can  also  be  used  to  apply  boundary  loads),  and  At 
is  the  time  step  size. 

Most  importantly,  peridynamics  implements  mater¬ 
ial  failure  by  separating  “bonds,”  i.e.  removing  nodes 
from  a  given  node’s  horizon  so  that  they  do  not  inter¬ 
act.  The  simplest  damage  model  in  the  literature  is  to 
remove  a  bond  based  on  a  bond  strain  measure,  the 
linearized  version  of  which  is  given  as 

 ^mn  ’  ^mn 
mn 

with  the  definitions  of  and  given  above.  A 
maximum  strain  can  be  given  as  a  constant,  giving  the 
failure  criterion 


^mn  ^  ‘^max, 


(9) 


where  ^'max  is  a  maximal  bond  strain,  typically  given  in 
terms  of  a  fracture  energy  as 


/47rGo 
V  9E8  ’ 


(10) 


where  Go  is  an  energy  release  rate  (Silling  and  Askari 
2005).  If  Eq.  (9)  is  met  for  a  given  bond  pair  {m,n], 
node  n  is  removed  from  node  m’s  family  (and  vice 
versa),  and  not  included  in  the  summation  of  Eq.  (6). 
After  a  bond  breaks,  its  supported  load  is  transferred  to 
other  bonds  in  the  nodes’  families,  encouraging  more 
bond  breakage  and  leading  to  progressive  crack  forma¬ 
tion. 

While  peridynamics  provides  a  description  of  the 
elastic  wave  equation  that  is  accurate  under  a  con¬ 
verging  horizon  and  its  discretized  form  gives  a  useful 
damage  model,  in  general  its  dispersion  characteris¬ 
tics  are  not  consistent  with  a  linear  elastic  material. 
Other  numerical  methods  can  provide  more  accurate 
wave  propagation  modeling  and  augment  peridynam¬ 
ics’  superior  handling  of  discontinuities. 


2.2  Finite  difference  method 

A  finite  difference  formulation  follows  standard  elas¬ 
ticity  theory  and  discretizes  its  differential  operators 
with  finite  difference  stencils.  The  linear  elastic  wave 
equation  can  be  stated  as 


,_„  =  V.,r+b, 


(11) 
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where  bold  type  with  an  overbar  represents  a  second- 
rank  tensor,  b  is  a  body  force,  the  stress,  a  is  given 
by 

a  =  C  :  €  =  C  :  ^  (Vu  -h  uV) ,  (12) 


two  overbars  represent  a  fourth-rank  tensor,  and  C  is 
a  fourth-rank  constitutive  tensor.  In  3D,  the  isotropic. 


linear  elastic  constitutive  tensor  is 

given  by 

2/x  -|-  X 

A. 

A 

0 

0 

0  " 

X 

2/x  “h  A 

A 

0 

0 

0 

C  = 

X 

A 

2/x  “h  A 

0 

0 

0 

0 

0 

0 

2/x 

0 

0 

0 

0 

0 

0 

2/x 

0 

0 

0 

0 

0 

0 

2/x  _ 

(13) 

where  the  Lame  parameters  A.  and 

/X  are 

k  = 


Ev 

(1  -h  v)  (1  -  2v)’ 


and 


(14) 


E 

^  2(l-hv)’ 


(15) 


Now,  the  above  set  of  equations  can  be  specialized 
with  the  restrictions  on  Poisson’s  ratio  inherent  in  peri- 
dynamics.  For  example,  for  2D  plane  strain  the  consti¬ 
tutive  tensor  becomes 


C  = 


2 

-E 

5 


"3  1  0" 
1  3  0 
002 


(16) 


resulting  in  the  coupled  equations 
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2  / 

f  9^ 

9^ 

5^1 

dxdy 

92 

2  / 

(  9^ 

92 

5^1 

3  ^  t/y  4-2 

9y^  ^  dxdy 

\ 


Ux-\--^Uy\+by. 


(17) 


The  remaining  two  formulations  follow  similarly,  for 
example,  the  plane  stress  formulation  is  identical  but 
with  the  factor  of  2/5  £"  replaced  with  3/8 

Often,  finite  difference  discretizations  for  elasticity 
will  use  two  or  four  staggered  grids,  solving  for  the  dis¬ 
placement  on  one  grid  and  the  stress  on  the  second,  or 
different  components  of  each  on  four  grids  (Virieux 
1986).  (This  approach  is  known  as  the  Yee  cell  in 
electromagnetics  (Yee  1966).)  Though  staggered  grid 
approaches  can  be  more  accurate,  here  we  use  a  sin¬ 
gle  grid  to  ensure  compatibility  with  the  peridynamics 


grid.  (In  an  extension  of  this  method,  state-based  peri¬ 
dynamics  could  be  implemented  in  a  dual,  staggered- 
grid  discretization  with  the  more  standard  finite  dif¬ 
ference  formulation.)  Rather,  the  second  order  spatial 
derivatives  will  be  discretized  with  a  standard  central 
difference  approximation  given  as: 


9^ 


/  (x  +  Ax,  y)  -  2/  {x,  y)  +  f  {x-  Ax,  y) 
Ax^ 


(18) 


dy 


f  {x,y  +  Ay)  -  If  (x,  y)  +  f  {x,y  -  Ay) 
Ay'^ 


(19) 


and 

92  f  {x  +  Ax,y  -\-  Ay) 

— — /  (x,  y)  K  - — — - 

dxoy  Ax  Ay 

_  f{x  +  Ax,y-  Ay) 

Ax  Ay 

fix  -  Ax,  y  -  Ay) 
Ax  Ay 

_  f  jx  -  Ax,y  +  Ay) 
AxAy 


(20) 


where  Ax  and  Ay  are  the  discretization  sizes  in  v  and 
y  respectively.  The  temporal  discretization  is  again  the 
velocity  Verlet  method  shown  in  Eq.  (7). 


2.3  Numerical  dispersion 

Next,  we  discuss  the  numerical  dispersion  inherent  in 
both  peridynamics  and  EDM,  as  any  interface  between 
the  two  methods  with  different  dispersion  characteris¬ 
tics  will  result  in  wave  refiections.  Previous  work  has 
shown  that  the  peridynamics  formulation  displays  dis¬ 
persion,  even  before  discretization  (Weckner  and  Abe- 
yaratne  2005)  and  it  is  well  known  that  FDM  suffers 
from  numerical  dispersion.  In  ID,  the  linearized  form 
of  peridynamics  behaves  as  a  continuous  sum  of  differ¬ 
ences,  much  like  discretized  finite  difference  methods. 
For  example,  specializing  Eqs.  2  and  6  to  ID  gives 
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which  can  be  discretized  as 
9  ^  Un-  Urn 


(22) 


A  central  difference  scheme  can  be  recovered  by  setting 
the  horizon  to  be  equal  to  the  spatial  discretization  size 
8  =  Ax,  resulting  in  Vmn  =  Ax,  Nm  =  2  (the  two 
nodes  adjacent  to  node  m),  c  =  EjS^,  and  giving 


Z 
/!=— 1 


E  U„-Um  U„-i  -2Um  +U„+i 


(23) 


A  standard  discretization  of  ID  peridynamics  with 
a  horizon  size  equal  to  the  spatial  discretization  size 
is  then  equivalent  to  a  central  difference  scheme,  thus 
motivating  an  augmented  method.  With  this  discretiza¬ 
tion,  the  best  that  can  be  done  is  to  match  the  dispersion 
of  a  finite  difference  technique,  though  it  will  be  clear 
that  including  a  larger  horizon  size  impacts  dispersion. 
First,  the  dispersion  relation  for  a  second  order  accurate 
in  space  and  time  FDM  can  be  given  by 

cos  (coAt)  =  y^cos  (PAx)  —  H-  1,  (24) 


where  y  =  vAt/Ax  is  the  Courant  number,  i;  = 
y/E/p  is  the  wave  speed,  (o  is  the  angular  frequency 
and  P  is  the  wavenumber. 

In  1 D,  the  discretized  peridynamic  operator  becomes 
a  sum  of  central  differences  in  all  of  the  nodes  within 
the  horizon.  The  addition  of  these  points  results  in  a 
different  dispersion  relation  for  each  value  of  the  hori¬ 
zon  size  8.  For  example,  with  a  horizon  size  of  2 Ax 
and  a  central  difference  discretization  in  time,  we  have 


cos  {(oAt) 


cos  (PAx)  -h  ^  cos  (2^21jc)1 


(25) 


and  for  a  horizon  size  of  3 Ax  we  have 
2y^ 

cos  {(oAt)  =  — ^  cos  (PAx) 

2  2 
-h  ^  COS  (ipAx)  -h  ^  cos  (3^^jc) 

-  + 1,  (26) 


As  an  example.  Fig.  1  gives  a  comparison  of  the  disper¬ 
sion  curves  of  FDM  and  peridynamics  with  two  horizon 
sizes,  8  =  2Ajc  and  8  =  3 Ax  for  Ax  =  1  mm,  At  = 
10  ns,  E  =  65  GPa,  and  p  =  2,235  kg/m^. 


Fig.  1  Dispersion  curves  for  FDM  and  discretized  peridynamics 
with  two  horizon  sizes 


Again,  as  we  wish  to  apply  FDM  to  some  areas  of 
the  computational  domain,  and  peridynamics  to  others, 
these  interfaces  will  exhibit  frequency-dependent  wave 
reflections  because  of  the  differing  dispersion  relations. 
The  magnitude  of  these  reflections  will  be  investigated 
in  Sect.  2.4.1. 


2.4  Augmented  method 

The  augmented  method  essentially  implements  a 
straight-forward  FD  discretization  while  simultane¬ 
ously  computing  peridynamics-based  bond  strains  as  a 
fracture  criterion.  Unless  fracture  occurs,  peridynamics 
is  not  used  to  compute  internal  forces,  it  is  only  used  as 
a  fracture  criterion.  After  bond  breakage  occurs,  peri¬ 
dynamics  is  then  used  to  compute  the  internal  forces. 

2,4.1  Wave  reflection  at  PD/FDM  interfaces 

As  mentioned  in  Sect.  2.3  FDM  exhibits  numerical  dis¬ 
persion  and  peridynamics  has  inherent  dispersion  prop¬ 
erties,  which  can  be  exacerbated  when  discretized.  Our 
computational  grid  will  contain  areas  that  use  both 
FDM  and  peridynamics  to  compute  wave  propagation, 
and  so  any  interface  between  the  two  will  result  in 
reflections.  The  reflection  coefficients  can  be  computed 
based  on  the  dispersion  relations  given  in  Eqs.(24), 
(25),  and  (26)  using  the  numerical  phase  velocity.  Here, 
we  will  simply  give  the  relation  in  ID,  which  would  be 
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representative  of  a  normally  incident  plane  wave  in  an 
infinite  medium.  The  reflection  coefficient  between  two 
materials  with  different  impedances  is  given  by 


Z2-Z, 
Z2  +  z,  ’ 


(27) 


where  Z\  and  Z2  are  impedances  of  the  two  media.  In 
1 D,  the  acoustic  impedance  is  given  by 


Z  =  =  pv. 


(28) 


If  we  use  the  phase  velocity  Vp  =  a)/P  and  the  dis¬ 
persion  relations  from  Eqs.  24, 25,  and  26,  we  can  sub¬ 
stitute  a  numerical  impedance  into  Eq.(27),  resulting 
in 

(29) 

(Opm  +  ^FDM 

where 

=  cos'*  l^y^cos(^^x)  ~  -b  l]| 

(30) 


and 


(opiy2  —  cos 


At 


-I 


[«■ 


cos  iPAx)  -h  -  cos  up  Ax) 


]|- 


(31) 


Similar  results  follow  for  other  combinations  of  disper¬ 
sion  relations,  three  of  which  are  plotted  in  Fig.  2  using 
the  same  parameters  as  those  in  Fig.  1 .  As  expected,  the 
wave  reflection  between  computational  grids  becomes 
worse  with  increasing  wavenumber  and  horizon  size; 


Fig.  2  Reflection  coefficients  for  normal  incidence  between 
FDM  and  discretized  peridynamic  computational  grids 


however,  this  is  only  an  illustrative  example  as  in  2D 
and  3D  the  reflection  coefficient  will  differ  for  angle  of 
incidence  and  wave  type  (longitudinal  vs.  shear).  We 
do  of  course  expect  similar  behavior  in  general:  Higher 
reflection  coefficients  for  higher  wavenumbers  (or  fre¬ 
quencies).  Finally,  while  larger  horizon  sizes  result  in 
larger  reflections  between  the  FD/peridynamics  inter¬ 
faces,  it  has  been  observed  that  larger  horizons  yield 
more  accurate  crack  features  (Ha  and  Bobaru  201 1). 
Ultimately,  a  large  enough  horizon  should  be  used  to 
model  fracture,  though  continually  increasing  the  hori¬ 
zon  will  result  in  less  accurate  behavior  at  computa¬ 
tional  interfaces.  Consequently,  in  the  dynamic  fracture 
results  (Sect.  3.4)  a  horizon  of  32\jc  will  be  used. 

2.4.2  Methodology 

The  augmented  method  proceeds  as  follows:  The  set 
of  nodes  is  split  into  two  sets,  nodes  at  which  inter¬ 
nal  forces  are  computed  using  finite  differences  Apo 
and  those  at  which  internal  forces  are  computed  using 
peridynamics  Apo-  These  node  sets  do  not  overlap 
(Afd  n  ApD  =  0)  and  together  they  comprise  the 
set  of  all  nodes  A(Afd  U  Apd  =  N).  At  each  node, 
I’m  ^  both  the  finite  difference  approximation 

of  the  internal  force  (divergence  of  the  stress)  and 
the  peridynamic  bond  strains  are  computed,  restated 
here  as: 

"(32) 

where  and  Uy  are  functions  of  r  and  t,  and  with 
the  appropriate  finite  differences  applied  to  each  dif¬ 
ferential  operator.  Simultaneously,  at  each  node  in  Apo 
the  peridynamic  bond  strain  is  computed  as  given  in 
Eq.  (8),  though  it  is  not  included  in  the  internal  force 
estimation  until  the  failure  criterion  is  met  (see  Eq.  9). 
After  some  bond  strain  Smn  exceeds  the  failure  criterion, 
the  internal  force  estimation  switches  to  peridynamics, 
again  given  as 

F™  =  Z  V„„,  (33) 

thus  moving  the  node  from  the  set  Apo  to  the  set  Apo- 
In  this  way,  peridynamics  is  used  similar  to  cohe¬ 
sive  zones  in  finite  element  methods:  After  some  field 
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measure  exceeds  a  critical  value,  a  new  constitutive 
response  is  used  to  model  the  internal  force  compu¬ 
tation  giving  progressive  failure/fracture  modeling.  In 
this  way,  the  finite  difference  force  computation  can 
more  accurately  model  the  wave  propagation  in  the 
bulk,  and  peridynamics  can  be  used  in  damaged  areas. 
This  approach  is  summarized  in  Algorithm  1 . 


Data:  E,  p,  8,  Ax,  Ay,  ^max 

Result:  u 

compute  integration  weights  at  each  node  for  peridynamic 

horizon,  8 ; 

compute  finite  difference  coefficients; 

for  each  time  step  do 

for  each  node  m  in  Afd  do 

compute  force  using  finite  differences; 
compute  bond  strains,  Smn ; 
if  any  Smn  >  ^^max  then 
break  bond  {m,n}\ 
switch  node  m  to  peridynamics  force 
computation; 
end 
end 

for  each  node  m  in  Apd  do 

compute  force  using  peridynamics; 
compute  bond  strains,  Smn ; 
if  any  Smn  >  ^^max  then 
break  bond  {m,  n}; 
end 
end 

velocity  Verlet  temporal  integration  to  update  u; 
apply  boundary  conditions  to  u; 

end 

Algorithm  1:  FDM  augmented  peridynamics  method 


2.4.3  Boundary  conditions 

Boundary  conditions  can  be  handled  using  either 
method,  peridynamics  or  finite  differences.  Dirich- 
let  boundary  conditions  are  straightforward  in  either 
method.  A  traction-free,  Neumann-like  boundary  con¬ 
dition  in  peridynamics  is  straightforward,  and  is  imple¬ 
mented  by  excluding  bonds  in  a  horizon.  In  other 
words,  a  traction-free  boundary  conditions  is  set  up 
when  bonds  are  broken  along  a  fracture.  In  the  finite 
difference  domain,  the  Neumann  condition  can  be 
set  by  expanding  the  traction  as  a  finite  difference 
of  boundary  and  internal  nodes,  and  solving  for  the 
unknown  boundary  displacements  in  terms  of  the 
known  internal  displacements.  For  example,  on  a  rec¬ 
tangular  domain  aligned  with  the  Cartesian  axis,  the 
traction-free  boundary  condition  on  the  top  edge  can  be 
stated  as 


y  a  =  0, 


(34) 


giving  the  coupled  equations 


9  9^ 

l/x  Uy  -  0, 

dy  dx  ^ 


d  d 

+3— My 

dx  dy 


0, 


(35) 


where  the  factor  of  3  in  the  second  equation  is  a  result 
of  the  constitutive  relation  of  Eq.  (16).  Applying  finite 
differences  to  the  above,  the  unknown  displacements  on 
the  boundary  can  be  written  in  terms  of  known  internal 
displacements  as 


Ayuy  (x  -h  Ax,  y)  —  Ayuy  (x  —  Ax,  y) 

-\-2AxUx  (x,  y)  =  2AxUx  (x,  y  —  Ay) , 

Ayux  (v  -h  Ax,  y)  —  Ayu^  (v  —  Ax,  y) 

-h  6AxUy  (x,  y)  =  6AxUy  (x,  y  —  Ay) ,  (36) 

where  a  central  difference  has  been  used  for  the  v  deriv¬ 
ative  and  a  backward  difference  has  been  used  for  the 
y  derivative.  A  sparse  matrix  can  be  formed  to  solve 
for  the  displacement  at  all  boundary  nodes  after  the 
displacement  update  in  Eq.  (7)  is  used  to  compute  the 
internal  displacements.  Of  course,  in  Eq.  (35),  the  zeros 
on  the  right  hand  sides  can  be  replaced  with  a  loading 
function,  taking  care  to  include  Young’s  modulus  and 
the  correct  scaling  factor  for  either  plane  strain  or  plane 
stress. 

Unfortunately,  it  has  been  found  that  any  inter¬ 
face  between  the  finite  difference  method  and  peri¬ 
dynamics  on  a  boundary  can  lead  to  high  enough 
stress  concentrations  that  spurious  fracture  can  form. 
For  now,  the  boundary  of  the  material  will  be  set 
to  not  accumulate  damage  other  than  pre-cracked 
regions. 


3  Numerical  results 

In  this  section,  several  numerical  results  will  be  pre¬ 
sented.  The  first  will  focus  on  comparing  the  dispersion 
effects  of  the  finite  difference  discretization  and  dis¬ 
cretized  peridynamics.  Next,  reflections  at  the  interface 
between  the  two  computational  regions  will  be  demon¬ 
strated.  In  Sect.  3.3,  the  stress  intensity  factor  of  a  semi¬ 
infinite  line  crack  in  an  infinite  medium  will  be  com¬ 
pared  in  both  the  hybrid  and  peridynamic  approaches. 
Finally,  a  dynamic  fracture  problem  will  be  illustrated 
in  Sect.  3.4. 


^  Springer 


46 


R.  A.  Wildman,  G.  A.  Gazonas 


3.1  Dispersion  comparison 

In  this  subsection,  the  dispersion  effects  of  a  finite 
difference  method  and  discretized  peridynamics  with 
two  horizons  {2Ax  and  3Z\jc)  in  a  2D  plane  strain 
setting  are  compared.  The  constitutive  parameters  are 
given  by  Young’s  modulus  £  =  65  GPa  and  density 
p  =  2,235  kg/m^  and  the  discretization  parameters  are 
given  by  jc-step  of  Ax  =  1  mm,  y-step  of  Ay  =  1  mm, 
and  time  step  size  of  At  =  10  ns.  Poisson’s  ratio  is 
fixed  at  V  =  1/4  as  discussed  above.  The  computa¬ 
tional  region  is  a  200mm-by-200mm^  (spanning  the 
region  0  <  jc,  y  <  200  mm),  chosen  to  be  large  enough 
so  boundary  effects  do  not  influence  the  results.  A 
plane  wave  with  a  Gaussian  profile  is  used  as  an  initial 
condition — so  as  to  compare  with  an  exact  solution — 
given  by 

U  (r,  t  =  0)=  exp  [-7.5  x  lO'*  (r  •  x  -  O.l)^]  x. 

(37) 

The  simulation  is  run  for  10p.s  and  Ux  is  observed  at 
the  points  r  =  (70,  100)  mm  and  r  =  (80,  100)  mm. 

First,  Fig.  3  shows  the  jc-directed  displacement  Ux 
vs.  time  at  the  two  points  listed  above.  The  displace¬ 
ment  measured  at  80  mm  is  shown  as  the  solid  lines 
and  the  displacement  at  70  mm  is  shown  as  the  dashed 
lines.  An  “exact”  solution  (determined  assuming  a  1 D 
problem,  though  the  numerical  solutions  were  com¬ 
puted  in  2D)  is  shown  as  the  red  curves,  the  finite  dif¬ 
ference  approximated  solution  is  shown  as  the  black 


Fig.  3  The  x-directed  displacement  versus  time  step  measured 
at  two  points,  demonstrating  dispersion  effects 


Frequency  (Hz) 


Fig.  4  Fourier  transform  of  Fig.  3 


Fig.  5  Absolute  error  in  the  Fourier  transform  of  Fig.  3  measured 
against  the  exact  solution 

curves,  and  two  peridynamics  approximated  solutions 
are  shown  in  green  and  blue,  using  horizons  of  3  Ax  and 
2 Ax  respectively.  As  can  be  seen,  peridynamics  with  a 
higher  horizon  displays  errors  characteristic  of  disper¬ 
sion  effects,  namely  a  lowered  amplitude  and  oscilla¬ 
tions.  Next,  the  Fourier  transform  of  the  Jc-directed  dis¬ 
placement  was  computed  and  is  shown  in  Fig.  4,  along 
with  the  absolute  error  in  each  solution  compared  with 
the  “exact”  solution  (Fig.  5).  These  figures  demonstrate 
that,  compared  to  the  finite  difference  approximation, 
higher  frequencies  in  the  peridynamics  methods  dis¬ 
play  larger  errors,  consistent  with  dispersion  errors. 
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3.2  Interface  reflections 

Wave  reflections  between  computational  regions  are 
examined  here.  This  was  discussed  in  an  analytical 
form  for  ID  in  Sect. 2.4.1,  though  here  we  perform 
numerical  tests  of  the  2D  algorithms.  In  the  first  test, 
the  parameters  of  the  problem  are  identical  to  those 
of  the  previous  subsection,  though  now  the  left  half  of 
the  domain  will  be  computed  with  the  finite  difference 
method  and  the  right  half  will  use  peridynamics  with  a 
horizon  size  of  either  3  Ax  or  2 Ax.  The  initial  condi¬ 
tion  is  again  a  Gaussian  pulse,  now  centered  at  the  line 
jc  =  75  mm  and  the  jc -directed  displacement  is  now 
measured  at  the  point  r  =  (90,  1(X))  mm. 

As  shown  in  Fig.  6,  the  reflected  wave  can  be  exam¬ 
ined  by  also  computing  the  solution  to  this  problem 
with  a  fully  finite  difference  domain,  and  subtracting 
the  results,  thus  leaving  only  the  reflected  wave  in  the 
left  half  domain.  In  Fig.  6,  the  blue  curve  is  the  result 
using  peridynamics  with  a  horizon  of  3  Ax  in  the  right- 
half  and  the  red  curve  shows  the  simulation  with  a 
horizon  of  2  Ax.  Given  that  the  magnitude  of  the  inci¬ 
dent  pulse  was  0.5  m,  the  reflections  at  the  interface 
can  be  significant  for  the  3  Ax  horizon  size.  By  max¬ 
imum  magnitude,  the  reflected  wave  is  approximately 
12%;  however,  peridynamics  will  only  ever  be  used 
in  regions  that  have  experienced  fracture  or  failure,  so 
some  amount  of  wave  reflection  would  be  expected 
regardless. 

Given  the  above,  it  may  however  also  be  possible 
that  the  anomalies  seen  in  Fig.  6  may  be  due  to  the 


Fig.  6  Absolute  error  between  a  fully  FD  solution  and  a  solution 
with  half  peridynamics,  half  FD 


Fig.  7  Absolute  error  between  a  fully  FD  solution  and  a  solution 
with  a  vertical  line  of  peridynamic  nodes 

dispersion  properties  of  the  peridynamic  region.  To  test 
this,  we  again  perform  the  above  simulation  in  a  FD 
domain  with  peridynamics  “window”  consisting  of  a 
single  vertical  line  of  nodes  at  the  line  jc  =  1(X)  mm. 
The  smaller  window  will  reduce  the  dispersion  effects 
of  a  wave  propagating  through  a  peridynamic  region, 
but  still  demonstrate  the  effects  of  interface  reflection. 
Figure  7  gives  the  results  of  this  test,  with  the  blue  curve 
showing  peridynamics  with  a  horizon  of  3  Ax  and  the 
red  curve  having  a  horizon  of  2 Ax. 

For  a  fully  fractured  region,  a  free  surface  is  formed 
and  \00%  of  a  wave  is  reflected,  so  as  a  final  test, 
we  introduce  a  free  surface  with  boundary  conditions 
implemented  using  finite  differences  to  compare  a  frac¬ 
ture  surrounded  by  peridynamic  nodes.  The  same  set  up 
is  used  as  described  above,  though  here  the  free  surface 
is  placed  at  jc  =  1(X)  mm.  The  initial  pulse  is  centered 
at  the  line  jc  =  70  mm  and  the  result  shown  in  Fig.  8  is 
measured  at  the  point  r  =  (90,  100)  mm.  Again,  Fig.  8 
shows  the  difference  between  the  fully  FD  simulation 
and  the  simulation  using  a  peridynamic  boundary  with 
a  thickness  of  one  horizon. 

3.3  Stress  intensity  comparison 

Next,  the  augmented  method  and  peridynamics  are 
compared  using  the  dynamic  mode  I  stress  intensity 
factor  at  a  semi-infinite,  stationary  crack  in  an  infi¬ 
nite  medium.  This  problem  has  an  analytical  solution 
for  the  stress  intensity  factor  for  both  stationary  and 
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Fig.  8  Absolute  error  between  a  fully  FD  solution  and  a  solution 
with  a  peridynamic  border 


Fig.  9  Normalized  stress  intensity  factor  comparison  between 
hybrid  and  peridynamics 


moving  cracks  (Liu  et  al.  2011),  though  here  we  focus 
on  the  stationary  case.  As  given  by  Liu  et  al.  (201 1), 
the  stress  intensity  factor  for  this  geometry  can  be 
stated  as 


Ki  it)  =  H{t-  fc) 


2cr*  lv\(t  —  tc)  (1  —  2v) 

1  -  vV  n 


(38) 


where  //  (•)  is  the  Heaviside  step  function,  ui  is  the 
longitudinal  wave  speed,  and  a*  is  the  applied  stress. 
To  simulate  an  infinite  problem,  a  large  region  is  used 
with  a  horizontal  pre-crack  through  the  center  and  the 
simulation  is  stopped  before  the  reflected  wave  reaches 
the  crack  tip.  The  material  properties  used  give  a  lon¬ 
gitudinal  wave  speed  of  v\  =  7,332  m/s  with  Pois¬ 
son’s  ratio  1/4.  Finally,  the  applied  tensile  stress  is 
40  MPa. 

To  estimate  the  stress  intensity  factor  in  the  sim¬ 
ulation,  the  method  outlined  by  Chen  (1975)  will  be 
used.  To  summarize,  the  stress  at  the  crack  tip  can  be 
expressed 


^,*(0/,  .  ^  .  3^\  9 


(39) 


where  K*  (f)  is  a  low-order  estimate  of  the  stress  inten¬ 
sity,  r  is  the  distance  to  the  crack  tip,  and  9  is  the  angle. 
The  stress  intensity  factor  Ki  (t)  can  then  be  estimated 
by  choosing  nodes  along  a  line  with  constant  angle  with 
respect  to  the  crack  tip,  computing  the  stress  at  those 
points  using  finite  differences,  and  constructing  a  least 
squares  approximation  of  the  stress  near  the  crack  tip. 


More  precisely,  we  choose  nodes  along  the  hne  ^  =  0 
(typically  five  nodes,  starting  with  the  fourth  closest  to 
the  crack  tip)  and  use  the  higher  order  approximation 
of  the  stress  intensity 

K*  ^Ki+a  (0)  r'/^  +  ^  (0)  r  +  y  (9)  (40) 

and,  according  to  Chen  (1975),  the  coefficients  of  non¬ 
integer  powers  of  r,  a  (9)  and  y  (9),  are  zero  for  9=0. 
A  least  squares  approximation  is  then  used  to  estimate 
Ki(t)  (Chen  1975).  Figure  9  shows  the  normahzed 
stress  intensity  factor  (normalized  by  Oy/na  where 
a  =  150  mm  is  the  crack  length)  with  the  exact  solu¬ 
tion  and  both  peridynamics  and  the  hybridized  method. 
The  oscillations  seen  in  peridynamics  (most  likely  due 
to  wave  dispersion)  are  reduced  in  the  hybrid  method. 

3.4  Dynamic  fracture 

Finally,  a  dynamic  crack  propagation  problem  is  simu¬ 
lated.  The  problem  setup  is  that  described  in  (Kalthoff 
1974):  A  150  mm  horizontal  crack  in  a  300mm-by- 
100mm-by-9mm  glass  plate  is  propagated  with  a  ten¬ 
sile  load  of  3,040  N.  No  material  properties  are  given 
in  (Kalthoff  1974),  so  we  use  properties  given  in  (Ha 
and  Bobaru  2010)  for  soda-lime  glass:  Young’s  mod¬ 
ulus  of  72GPa,  density  of  2,440  kg/m^,  and  Poisson 
ratio  of  i/3  (restricted  for  plane  stress).  We  also  have 
chosen  a  fracture  energy  value  of  3  J/m^,  which  gives 
a  crack  pattern  that  most  closely  matches  the  exp>er- 
iment  in  (Kalthoff  1974).  While  this  fracture  energy 
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Fig.  1 1  Magnification  of  the  crack  branching  area  of  Fig.  10 

appears  somewhat  low,  it  is  close  to  the  value  of  approx¬ 
imately  8  J/m^  found  in  the  literature  for  crack  initiation 
(Doll  1975).  Finally,  in  the  peridynamics  simulations, 
the  load  is  applied  as  a  body  force  as  described  in  Ha 
and  Bobaru  (2009),  and  as  a  traction  boundary  condi¬ 
tion  in  the  hybrid  method,  with  a  traction  of  1 . 13  MPa. 

First,  discretization  sizes  of  Ax  =  1  mm.  Ay  = 
1  mm,  and  At  =  10  ns  were  used  in  both  methods  and 
the  simulations  were  run  for  100  |xs  ( 10,000  time  steps) 
with  a  horizon  size  of  <5  =  3 Ax.  The  final  crack  pat¬ 
terns  are  shown  in  Fig.  10,  where  a  crack  is  taken  to 
be  any  node  with  damage  >0.4  (as  measured  by  the 
ratio  of  the  number  of  broken  bonds  to  the  number  of 
initial  bonds)  and  the  hybrid  result  is  shown  in  red  and 
the  peridynamics  is  shown  in  blue.  Figure  1 1  shows  a 
magnified  view  of  the  crack  branch:  As  can  be  seen, 
both  methods  branch  at  approximately  the  same  loca¬ 
tion,  though  they  differ  somewhat  in  their  subsequent 
crack  paths. 

The  same  problem  was  repeated,  though  now  with 
discretization  sizes  of  Ajc=0.5mm  and  Ay =0.5  mm, 
and  with  a  horizon  again  of  5  =  3Ajc.  Figure  12  shows 
the  results  of  the  simulations,  again  with  the  hybrid 
method  in  red  and  the  peridynamics  in  blue.  Here,  the 
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Fig.  12  Magnification  of  the  crack  branching  area  with  a  finer 
discretization 


peridynamic  method  alone  appears  to  have  an  earlier 
crack  initiation  time.  This  is  investigated  further  by 
tracking  the  crack  tip  position  and  speed  versus  time. 

The  crack  speeds  were  estimated  (for  the  finer  dis¬ 
cretization)  by  tracking  the  damage  as  it  propagated 
across  the  domain  and  are  shown  in  Fig.  13.  The  loca¬ 
tion  of  the  crack  was  recorded  at  every  5  |xs  (500  time 
steps)  (shown  in  Fig.  14)  and  the  crack  speed  was  esti¬ 
mated  by  computing  the  central  difference  of  the  result. 
The  hybrid  method  has  a  crack  initiation  time  slightly 
later  than  the  peridynamics  results  and  a  crack  speed 
slightly  slower  than  that  of  the  peridynamics.  These  dif¬ 
ferences  are  most  likely  due  to  the  lower  stress  intensity 
factor  observed  (see  Fig.  9)  at  the  crack  tip  of  the  hybrid 
method.  This  can  be  confirmed  by  comparing  the  end 
points  of  the  cracks  in  Fig.  12.  Two  features  can  be  seen 
in  the  crack  speed  data  in  the  form  of  the  two  dips  in 
crack  speed  around  35  and  65  p,s.  The  first  drop  in  speed 
occurs  as  the  crack  branches  and  the  second  when  the 
branch  angles  change  slightly  at  225  mm. 

It  is  also  illustrative  to  compare  both  methods 
between  the  two  discretizations  given  above.  To  that 
end.  Fig.  15  shows  the  peridynamic  method  at  both 
of  the  above  discretizations.  It  can  be  seen  in  the  fig¬ 
ure  that  the  crack  initiation  occurs  slightly  earlier  in 
the  finer  discretization  (shown  in  red).  In  contrast,  the 
hybrid  method  shows  a  more  stable  fracture  path,  with 
both  following  a  very  similar  path  and  initiation  time 
(Fig.  16).  In  both  figures,  the  coarse  discretization  is 
plotted  in  blue  and  the  fine  discretization  is  shown  in 
red.  These  differences  are  most  likely  caused  by  the 
oscillations  seen  in  the  stress  intensity  factor,  due  to  the 
higher  dispersion  in  the  peridynamic  method  alone. 
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Fig.  13  Comparison  of  crack  speeds  of  the  hybrid  method  {red) 
and  peridynamics  {blue) 


Time  (|is) 

Fig.  14  Comparison  of  crack  tip  locations  of  the  hybrid  method 
{red)  and  peridynsimics  {blue) 


Finally,  the  dynamic  fracture  problem  was  repeated 
with  both  discretizations  used  above,  but  with  a  higher 
load  of  6080  N.  In  this  case,  a  more  complicated  crack 
path  is  observed,  with  a  larger  difference  between  the 
methods.  Figure  17  shows  the  result  using  discretiza¬ 
tion  sizes  of  Ax  =  1  mm,  Ay  =  1  mm.  A  major 
difference  between  the  two  methods  is  the  shallower 
crack  branch  angle  at  the  branch  onset  and  the  addi¬ 
tional  branch  close  to  the  center  in  the  hybrid  method. 
Otherwise,  they  both  feature  multiple  branches,  con¬ 
sistent  with  crack  patterns  seen  at  higher  loads  (Bow¬ 
den  et  al.  1967;  Ha  and  Bobaru  2011).  Figure  18  shows 
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Fig.  15  Comparison  of  crack  paths  at  different  discretizations 
(coarse  in  blue  and  fine  in  red)  using  peridynamics 
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Fig.  16  Comparison  of  crack  paths  at  different  discretizations 
(coarse  in  blue  and  fine  in  red)  using  the  hybrid  method 
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Fig.  1 7  Comparison  of  crack  paths  at  a  higher  load  and  a  coarse 
discretization 

the  result  using  discretization  sizes  of  Ax  =  0.5  mm. 
Ay  =  0.5  mm.  At  this  discretization,  the  peridynam¬ 
ics  solution  features  an  additional  branch  near  the  far 
edge  and  two  arrested  cracks  not  seen  at  the  coarse 
discretization  of  the  peridynamics  method. 
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Fig.  18  Comparison  of  crack  paths  at  a  higher  load  and  a  fine 
discretization 

4  Conclusions 

A  hybrid  finite  difference/peridynamics  method  was 
presented  to  improve  the  dispersion  properties  of  stand¬ 
alone  peridynamics.  The  hybrid  method  essentially 
switches  the  spatial  internal  force  computation  between 
finite  differences  (divergence  of  the  stress)  and  peri¬ 
dynamics  (integral  of  a  micro-force  function)  after  a 
damage  criterion  is  met.  Linearized  peridynamics  can 
be  shown  to  model  linear  elastic  materials  with  specific 
Poisson  ratios,  and  this  fact  is  used  to  develop  a  match¬ 
ing  finite  difference  scheme.  Analysis  was  presented 
showing  the  dispersion  properties  of  discretized  peri¬ 
dynamics  in  ID  in  comparison  to  a  straight-forward 
finite  difference  discretization.  Numerical  results  were 
also  shown  to  demonstrate  the  dispersion  properties  of 
the  two  methods  in  2D.  Other  results  show  that  the 
numerically-computed  stress  intensity  factor  for  peri¬ 
dynamics  can  oscillate,  which  may  lead  to  differing 
crack  paths. 

This  method  has  a  number  of  possible  extensions, 
such  as  3D,  state-based  peridynamics,  and  non-linear 
materials.  Additional  work  may  be  done  as  well  to  ease 
the  transition  region  between  peridynamics  and  finite 
differences,  as  reflections  can  occur  due  to  the  differing 
material  representations. 
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